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Up to now, the effects of having heterogeneous networks of contacts have been studied mostly for diseases 
which are not persistent in time, i.e., for diseases where the infectious period can be considered very small 
compared to the lifetime of an individual. Moreover, all these previous results have been obtained for closed 
populations, where the number of individuals does not change during the whole duration of the epidemics. Here, 
' we go one step further and analyze, both analytically and numerically, a radically different kind of diseases: 

, those that are persistent and can last for an individual's lifetime. To be more specific, we particularize to the 

case of Tuberculosis' (TB) infection dynamics, where the infection remains latent for a period of time before 
' showing up and spreading to other individuals. We introduce an epidemiological model for TB-like persistent 

infections taking into account the heterogeneity inherent to the population structure. This sort of dynamics 
introduces new analytical and numerical challenges that we are able to sort out. Our results show that also for 
' persistent diseases the epidemic threshold depends on the ratio of the first two moments of the degree distribution 

so that it goes to zero in a class of scale-free networks when the system approaches the thermodynamic limit. 

PACS numbers: 87.23.Ge, 89.20.-a, 89.75.Fb 
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^ ; I. INTRODUCTION 

\ Disease spreading has been the subject of intense research since long time ago IH-Ilt]. Our current knowledge comprises 
c/2 ■ mathematical models that have allowed to better understand how an epidemic spreads and to design more efficient immunization 
and vaccination policies IH-d!]. These models have gained in complexity in recent years capitalizing on data collections which 
have provided information on the local and global patterns of relationships in the population |4j-t6[|. In particular, with the advent 
of modern computational resources and tracking systems, it is now feasible to contact-trace the way the epidemic spreads or at 
least to predict the paths that a given pathogen might follow. In this way, some of the assumptions at the basis of the theoretical 
models that were difficult to test -the backbone through which the diseases are transmitted- are now more accurately incorporated 
into epidemiological models iTl- jTlll . 

Strikingly, the systems on top of which diseases spread show common nontrivial topological and statistical properties lfl2l[l3ll . 
_ A large number of networks of contacts in real-world social, biological and technological systems have been found to be best 

■ described by the so-called scale-free (SF) networks. In SF networks, the number of contacts or connections of a node with other 
nodes in the system, the degree (or connectivity) k, follows a power law distribution, P{k) ^ k^'^ . Recent studies have shown 
that the SF topology has a great impact on the dynamics and function of the system under study lfl2ljl5ll . The reason is that, 

■ at variance with homogeneous or regular networks, SF architectures are a limiting case of heterogeneity where the connectivity 
^-H . fluctuations diverges if 2 < 7 < 3 as the system size tends to infinity (the thermodynamic limit). This means that there are 

ILJ ' nodes in the network which has an eventually unbounded number of connections compared to the average degree. Examples of 
such networks include the Internet lfl6l[l7ll . the world- wide-web (WWW) ifisll . food-webs, and metabolic or protein networks 

5— i ] In the context of disease spreading, SF networks lead to a vanishing epidemic threshold in the limit of infinite population when 
^ ■ 7 < 3 ll20M23ll . This is because the ratio {k) / (fc^) determines the epidemic threshold above which the outbreak occurs. When 
2 < 7 < 3, (fc) is finite while (fc^) goes to infinity, that is, the transmission probability required for the infection to spread goes to 
zero. Moreover, the previous result holds both for the Susceptible-Infected-Susceptible (SIS) and Susceptible-Infected-Removed 
(SIR) epidemiological models ll20l - l22ll . 

In this paper, we will deal with a different kind of diseases -those that are persistent in time and shows a latent period that can 
be as large as an individual's lifetime. Our first aim is to enlarge the epidemiological framework for complex networks reported 
previously for the SIS model |20il and proposed as well for the SIR model 12 U i221 by integrating the spreading dynamics of 
pers istent diseases within it. With this purpose, we consider a variation of the Susceptible-Exposed-Infected-Removed model 
II24I1 on complex heterogeneous networks. As we will see, this kind of dynamics introduces new challenges essentially different 
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to other successfully treated before. In particular, as the latent period is high enough, we shall work with an open system in 
which new individuals are born and others dead for causes not directly related to the spreading of the disease. We present novel 
analytical and numerical methods that allow us to obtain the epidemic threshold for this kind of disease in heterogeneous popu- 
lations. Moreover, we introduce a numerical method that is well-suited to deal with the kind of problems we face. Our results 
point out that also for persistent infections the virtually unbounded fluctuations of the degree distribution have an important 
enhancing effect on the epidemic spreading. 



II. THE MODEL 



To be more precise and without loss of generality, we will particularize our model to one case of persistent infection -the most 
threatening one which is tuberculosis (TB). TB is an old disease whose world-wide prevalence had been diminishing even before 
vaccination and prophylaxis strategies were firstly accomplished IzsI - IztIi . Its recent return in developing countries, mainly in 
Southeast Asia, have attracted renewed interest in it. The current world estimate of prevalence is about 33% while the number 
of deaths per year that it is causing reaches more than 3 million people ll28ll . Depending of the kind and the intensity of immune 
response that the host immune system performs after initial infection with M. tuberculosis bacillus, the individual can suffer 
latent infection, (in which the bacteria are under a growth-arrest regime and the individual neither suffer any clinical symptom 
nor becomes infectious) or active infection, where the host suffers clinical symptoms and can transmit the pathogen by air 
ll29l[30ll . Latently infected individuals can, generally after an immune-depression episode, reach the active phase. Estimating the 
probability of developing direct active infection after a contact, or alternatively, the lifetime's risk for a latent infected individual 
to evolve into the active phase, are not easy tasks. However, it is generally accepted that only 5-10% of the infections directly 
produce active TB ll29l [30ll . while the ranges concerning the estimation of typical "half-life" of latent state rounds about 500 
years (Ml- 

The spreading dynamics of TB like diseases has been studied in recent years. However, to the best of our knowledge, 
these works assume the homogeneous mixing hypothesis, that is, a perfectly homogeneous system in which all individuals are 
dynamically equivalent. As mentioned above, many of the systems on top of which diseases spread, are better described by 
scale-free connectivity patterns. Therefore, in what follows, our main objective will be to assemble a basic model fitted for 
tuberculosis spreading that firstly takes into account the heterogeneity in the distribution of the networks of contacts. We note 
that the increasingly alarming situation about TB epidemiology evidences the need to increase the effort in TB research in a 
global way. In the context of the study of its epidemiology, new models must be developed in order to gain predictive skills, 
incorporating the recent theoretical advances referring to disease spreading on complex heterogeneous substrates as well as 
meta-population approaches and new computational tools for numerical analysis and simulation. In this sense, ours is a first 
contribution that addresses one of the most important parameters in epidemiological description: the epidemic threshold. 

Let us then introduce our model. We consider that individuals in the population are compartmentalized into three groups: 
healthy -U{t)-, infected but not infectious -or latently infected L{t)- and sick individuals T{t) which are infected and are 
infectious as well. The transition between these subpopulations proceeds in such a way that a healthy individual acquires the 
bacteria through a contact with an infectious subject with probability A. In its turn, this newly infected individual may develop 
the disease directly with probability p. However, the most common case is the establishment of a dynamic equilibrium between 
the bacillus and the host's immune system, which allows the survival of the former inside the latter. When this happens, newly 
infected individuals become latently infected, because despite harboring the bacteria in blood, neither becomes sick nor is able 
to infect others. 

On the other hand, after a certain period of time (which may be several years) and usually following an episode of immuno- 
suppression, the balance between the bacterium and its host can be broken. In this case, the bacteria grow and the individual 
falls ill beginning to develop the clinical symptoms of the disease. In addition, if the infection attacks the lungs (pulmonary TB), 
the bacillus is present in the sputum, making the guest infectious. 

The dynamics of the disease, in a well mixed population, is then described by the following system of nonlinear differential 
equations: 

^ = (1 - P)m{t)§^^ - (m + r)m, (1) 

in which: 
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N{t) = U{t) + L{t) + T{t) represents the total population at time t, 
(3 is the number of contacts per time unit, 

A is the probability that the bacteria is transmitted to a new host after a contact with an infectious subject, 

h is the birth rate per capita and per unit time, 

jjL is the natural death rate per capita and unit time, 

jjLtb is the rate of disease-related deaths per capita and unit time, 

r is the transition frequency of latent infection (i.e., the probability that a latently infected individual becomes infectious), 
with the closure relationship: 

'^E^ = (h~p)N(t)-p.,,T{t). (2) 
at 

The model above is a variation of the archetypal SIR model to which a fourth class has been added (latency class L). This kind 
of model has been largely treated in the literature in its well-mixed version, and it is frequently referred as SEIR model. As a 
first step, in this work, we will identify the removed individuals mostly with dead ones, and therefore we do not consider the 
possibility of natural or medical recovery (this simplification is in part justified by the large latency period of infected individuals 
and the constant flow of newborns into the system). A more refined model would consist of introducing such eventual recovery 
fluxes in the model, as well as the possibility of further relapses (the so called endogenous reactivation). These phenomenologies 
might be important mainly for diseases (like TB) for which the only feasible treatment in many areas consists of supplying large 
series of antibiotics. Thinking on the tuberculosis case, further refinements, like the inclusion of varieties of less infectious 
extra-pulmonary diseases llsill . could also have consequences on the disease's dynamics. 



III. STRUCTURED POPULATIONS 



A. Dynamics 

The previous system of differential equations describes the dynamics of the epidemics in the well-mixed case. However, as 
argued above, the number of contacts of a given individual in a population can vary, which is reflected in an heterogeneous 
distribution of the number of contacts in the system. To account for this fact, we next consider a structured population described 
by a connectivity distribution P{k). The system of Eqs ([T]i has to be modified accordingly. Assuming that all individuals with 
the same number of contacts, i.e., belonging to the same connectivity class k, are dynamically equivalent, the new system of 
differential equations are formulated for each degree class. Therefore, for a structured population, we have that: 

Nk{t)^P{k)N{t), (3) 

with: 

Uk{t)+Lk{t)+n{t)=Nk{t). (4) 
Moreover, it is convenient to express the previous equations in terms of densities, also defined within each connectivity class: 

Uk{t) 



Uk{t) = 



Nuity 



m = (5) 



Nk{ty 

Tk{t) 



Nk{ty 

so that the following closure relation for any value of k is verified: 

Uk[t)+lk{t)+tk{t) = l y{k,t). (6) 
On the other hand, the probability O that any given link points to an infectious individual is given by: 
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which leads to the following set of equations that describes the dynamics within each connectivity class: 

dUk{t) 



dt 

dLkjt) 
dt 

dTkit) 



dt 

Finally, the number of individuals with connectivity k evolves according to: 

dNk{t) 



= bP{k)N - XkQ{t)Uk{t) - ^Uk{t), 

= (1 - p)\kQ{t)Ukit) - + r)Lu{t), (8) 
: p\k@{t)Uk{t) + rLk{t) + Htb)Tk{t). 



dt 



= {b- ti)Nk{t) ~ ^itbTk{t) Htbtk)Nk{t). (9) 



At this point, and building on the previous equation, it is important to point out a feature of the model: the influence of the 
infection dynamics on the connectivity distribution P{k). First, if we add the above equation for all k, we obtain that the total 
population evolves as: 



'^^^^^ {b f,)Nit) - ntb 



Y^^kit) = lb- ^i- ^itbJ2Pik)tk \ Nit). (10) 

fc \ k J 



dt 

However, if we substitute Nk{t) = P{k)N{t) directly into Eq. (|9]l and assume P{k) to be constant, we would arrive to: 

= P{k) (& - /i - fHbtk) N{t). 

The last expression is only compatible with Eq. ( fTol i under the unrealistic assumption that all connectivity classes have the 
same proportion of sick individuals. We must therefore assume that the distribution of connectivity is also a function of time: 
P(fc, t), and therefore: 

dNkjt) _ d[Pik,t)N{t)] _ ^^^^^^ dP{k,t) ^ p^j^^^^^dNjt) ^ ^jj^ 
dt dt dt ^ dt ' 

so, if we substitute Eq. (fTTT i into Eq. ^ we get: 

N{t)'^^%^ + P{k, t)^ = P{k) (fe - M - f^tbtk) N{t), 

expression from which, if we replace dN{t)/dt from Eq. (fTOl i. we get the temporal evolution of P{k, t) as: 

dP{k,t) 
Jt 

where 



-P{k,t)fitb[tk{t) - {tk){t)], (12) 



(ifc)W (13) 

k 

Reformulating the equations in terms of densities using the definitions of the densities given above, Eqs. (O become 

= b - uk{t)ib + xke{t) - ^Hbtk{t)), 

at 

= (1 - p)XkQ{t)uk{t) -{b + r)lk{t) + iitblk{t)tk{t), (14) 



dlkjt) 
dt 

dtk{t) 
dt 



: p\kQ{t)uk{t) + rlk{t) - {b + fitb)tk{t) + ^Hbtkitf 



5 



B. Evolution of the degree distribution 

At this point it is appropriate to point out one aspect that will hinder any numerical characterization of the epidemic threshold. 

Our main goal will be to calculate both numerically and analytically the critical value Ac beyond which the population presents 

Ik) 

an endemic proportion of sicks individuals. However, we expect Ac to be dependent on the ratio -^^y, which is in turn a function 
of the connectivity distribution P{k). The degree distribution, as previously shown, changes in time as the dynamics of infection 
progresses. As we should see, we can handle this time dependence analytically, but we should be forced to design a simulation 
method to account for the rate of births and deaths and the effects of these two processes on the degree distribution. 

The aforementioned features might lead to a situation in which the infection dynamics would modify the underlying structure 
of the network through which the disease is being spread. Therefore, Ac could also vary as one expects it to be intrinsically 
related to the first two moments of a seemingly time-dependent degree distribution. The reason why we consider the distribution 
of contacts per unit time as heterogeneous, even for the current airborne-transmitted disease is based on the observation that 
the number of contacts a person can have per unit of time is subjected to two sources of heterogeneity. Firstly, what we can 
call geo-demographic, macroscopic heterogeneity, in which the number of contacts depends on the population density in the 
region in which an individual inhabits. Secondly, at a more individual, microscopic level, the heterogeneity arises because the 
number of contacts depends, in a region of constant population density (i.e., a town or neighborhood in a city), on the daily 
activity pattern of the individual within that region. These two factors define, ultimately, the function P{k). Having that said, 
the assumption implicitly incorporated in the first equation of system (|8]l does not hold. Note that this equation implies that the 
connectivity of individuals is hereditary and therefore that the number of births within each k class equals the birth rate times 
the number of individuals within each k class, Nk = P{k, t)N . 

The above situation would be equivalent to assume that the dynamics of the disease being studied is the only one that influences 
the demographic structure of a population, which is not true since it is clear that there are countless cultural, economic and social 
factors that ultimately define the above two levels of heterogeneity. We therefore assume in what follows that the newborns of 
each generation are distributed among the k classes according to an invariant distribution function, which we further assume 
to be the initial degree distribution of the original network: P{k,to)- As we shall see, this assumption, besides being more 
plausible, has the advantage that makes the connectivity distribution to be roughly stable, and so will be the critical value Ac. 

So, we have the following reformulation of the system of differential equations (Ull: 

dUkit) ^ j,p(^^^^)^ _ \kQit)Ukit)-fiUkit), 



dt 

dLkjt) 
dt 

dTk{t) 



(1 - p)XkQ{t)Uk{t) + r)Lk{t), (15) 



: p\kQ{t)Uk{t) + rLkit) + Htb)Tk{t), 
dt 

with the definition of the number of individuals in each class of connectivity: 

Nkit) = N{t)Pik,t), (16) 

and inside each class: 

Uk{t) ^ Nkit)uk{t), 



Lk{t) = Nk{t)lk{t), (17) 



Tk{t) = Nk{t)tk{t). 
Now the total population within each connectivity class verifies: 

dNk{t) 



dt 



= bNP{k, to) - ^iNP{k, t) - nMt), (18) 



so that, if we add in k, the last modification has no effect on the variation of the total volume of the population. The temporal 
evolution of the degree distribution is now given as: 

^^^^ = b[Pik,to) - P{k,t)] - Pik,t)fitf, [tk{t) - {tk){t)] . (19) 
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Finally, writing the equations in terms of the densities we get 



dt ' ' V y -^v / y P(fc,t) 

+ ^itblk{t)tk{t), (20) 



C. Characterization of the equilibrium points. 

The previous set of differential equations tells us how the different densities of interest evolves within each connectivity class. 
Their corresponding macroscopic quantities are defined as 

{u){t)^Y.^P{k,t)uk{t), 

{l){t) = Y,kP{k,t)h{t), (21) 

{t){t) ^Y.^p{k,t)tu{t), 

where {u){t), {l){t) and {t){t) are the mean densities of healthy, latent and sick individuals, respectively. 

Let us now go one step further and characterize the equilibrium points. The magnitudes of interest are the average densities, 
so that an equilibrium point ((u)*, (/)*, (t)*) must verify by definition: 



d{u)* d{l)* d{t)* 



(0,0,0). 



dt dt dt 

We also impose a further constraint which is that the degree distribution of the network is stationary, that is: 

dP{ky 



dt 



Vfc. 



At this point one must ask whether macroscopic stability also implies stability within each connectivity class. The answer is 
yes, if we also demand stability of the degree distribution. Admittedly, if we equate expression (T9[ to zero and solve for the 
stationary P{k, t)* we get: 

Pik,tr^ 



b + fitb{tl- (t)*)' 



which shows that this value depends on the microscopic scale tk- Therefore, the stability of the degree distribution imposes a 
stationary condition on tk for all fc, which in its turns extends to the other densities l'^ and u^. Hence, we have: 



dul dll dtl 
dt ^ dt ^ dt 



(0,0,0) Vfc. 



The above condition is trivially satisfied for the solution {u*f, ,ll,tl) = (1,0, 0) Vfc, which leads to a degree distribution exactly as 
the initial distribution. We next analyze the stability of this solution, which shall allow us to characterize the epidemic threshold. 



D. Epidemic threshold 

As stated before, in this section we will study the stability of the solution {u*^, ll,tl) = (1, 0, 0) Vfc. At this point, as no latent 
or infected individuals are introduced in the network, the degree distribution does not change in time; so that P(fc)* = P{k, t) = 
P{k, to). This situation allows to work with the system of differential equations given by ( fT4b instead of working with the more 
general case given by the system ( l20l i. 
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1. Case p — 1 

For simplicity and to gain some preliminary insight into the problem, we first study the case p =\, which means that there is 
no latent phase (i.e., the latent subpopulation disappears, 1^ = Vfc). Using wj. + tj. = 1 we get, 

= h-Uk{h + Afce - /.itb) - /itfctfc, 
where we have omitted temporal dependences, as we will do from now on. Looking for the stationary solution, we have that the 

dt 

{ 1 



condition = impUes: 



from which the meaningful solution is the one with the negative sign. The previous expression is consistent with the meaning of 
u* since we recover the expected result u* ~ 1 when = 0. Moreover, if we calculate the derivative with respect to we get: 



dul{Q) _ Xk I ^ b + Xke~fitb 



<o, 



which guarantees that u* will always be less than unity and therefore is a real, valid solution. The study of the value of Q in the 
steady state help us to identify the epidemic threshold. We write: 



which, after substituting for its value, leads to: 



- :r^77Y y. kPik)^(bTxke~iIafTihhb. 

2ntb{k) ^ 

The graphical interpretation of the above equation indicates that the existence of an equilibrium point in which 9* > is 
equivalent to the existence of a point at which /(0) crosses the bisector of the first quadrant. Evaluating the second derivative 
of /(O) one gets: 

^V(Q) _ -2&A^ Y" P{k)k^ 

de^ ^ (fc) ^ [{b + XkQ - fitb)^ + Abfitb] ^ ' 

which ensures that the condition for the existence of such intersection is reduced to: 

f df{e) \ A(P) 1 _ 

V dQ /q^o (fc) b + ^tb 
from which the epidemic threshold is derived as: 

_{b+jMb){k)_ 
(P) • 

Note that apart from the factor {b + fitb), the previous result, formally coincides with the epidemic threshold of the SIR model. 

2. Case p 7^ 1 

This is a somewhat more involved case. For structured populations, the resolution of the system of differential equations 
(|20] | cannot be done explicitly. We next find the epidemic threshold for the case p ^ 1 using two approaches. On one hand, 
we study the time derivative of Q. On the other hand, we will also make use of the singularity of the Jacobian at the point 
{uk, Ik, tk) = (1, 0, 0) to argue that the expression for the critical threshold is given by: 

_ (fc) ir + b){ntb + b) 
^^>' " W) • ^^^^ 
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a. Time evolution of& in populations with a low number of sicks. A first approach to characterize the epidemic threshold 
in heterogeneous networks when p ^ 1 is to study the sign of the derivative of 9 at the onset of an epidemic outbreak. We 
consider an initially healthy population in which a small proportion of infectious individuals is introduced so that tk « 1 Vfc. 
The derivative of 6 is; 



dQ 
'dt 



dP(k) 
dt 



{k) 



Y.kPik)ktk 

{k) 



dP(fc) 
dt 



(k) 



which, after substitution of the values of the derivatives of P{k, t) and tk{t) leads to: 



/dB' 



Y.kP{k)kh 



pA0 



(6 + /itt,)e + /itt,e2. 



At this point we make two simplifications. The first and most easily justifiable is to neglect the term 8^. The second is related to 
the presence of 1^ in the above equation, that we have to transform in a dependency with respect to t^. Specifically, we assume 
to be sufficiently close to the stationary point (itfe, Z^, ife) = (1, 0, 0) as to be able to assume that the three derivatives vanish. In 
other words, and focusing our attention on latent and sick classes, we assume that: 

dlk 
dt 



-f = (1 - p)\kQuk -ib + r)lk + fitbhtk ^ 0, 



e~o 



from which: 



dt^ 
dt 



■ pXkQuk + rlk - (& + lJ.tb)tk + fJ-tbtl ~ 0, 



r+pb- iitbPtk 
which allows to express the derivative of Q as: 



{l-p)ib + fitb)tk-il-p)fMbtl {1 - p)ib + fitb) , , ^/,2n 

- tk + U(t,.), 



dt J e^o 



r [l - p){b + inb) 
r+pb 



r + pb 



{b + ntb)+ 



+pX 



Y^kPrnw 



(fc) 



In the limit Uk 
derived as: 



1 Vfc the third term within brackets is the ratio (fc^) / (fc), from which the epidemic threshold condition may be 



dB\ 



= e 



r {I - p){b + ptb) 
r + pb 



- ib + IHb) + P\ 



(fc2 



(fc) 



= 0, 



finally leading to the expected expression for the threshold: 

(fc) {r + b){^ltb + b) 



(fc2) pb -\- r 



(23) 



b. Analysis of the Jacobian While for well-mixed populations the condition of singularity of the Jacobian allows to get 
the epidemic threshold in a straightforward way, for heterogeneous populations the analysis of the Jacobian is a difficult task 
because is a function of each and every one of the i^'s. This translates into the need of calculating a determinant whose order 
is three times the number of connectivity classes. What we can reasonably do is to verify is the threshold condition is verified 
for systems in which there are two or three different connectivity classes. In the first case in which only two different classes of 
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FIG. 1 : Set of allowed transitions in the epidemic model within each connectivity class. 



connectivity exist, the Jacobian is just a quite distasteful 6x6 determinant that, after some cumbersome and lengthy algebra, can 
be reduced to the expression: 



{b + r){h + iJitb) + ^{pb + r) 



which equated to zero leads again to the previously obtained expression for the epidemic threshold. If we instead consider 
a population with three degree classes, the algebraic complexity of the problem largely increases as we now have to solve a 
determinant of size 9x9. However, we can proceed as before getting the following expression for the Jacobian: 



This leads us to the sensible conjecture that increasing the number of connectivity classes does not add new roots to the Jacobian, 
but it only would increase the degeneracy of the non-interesting solutions b = 0, b = —r and b = —fitb- 



E. Numerical simulations 



When designing numerical simulations to inspect the dynamics of the system under study, we have two difficulties not previ- 
ously addressed in the literature. These numerical issues with which we have to deal come from the fact that we have a system 
that is simultaneously open and structured. As a result of dealing with an open system, new individuals are being added to 
the population at a rate given by the birth rate. Additionally, these new individuals must enter the network of contacts with 
a predefined connectivity. While deciding how many nodes our new individuals connect to is not a problem, it certainly is to 
decide what are those nodes the newcomers will be linked to, as this will impact the degree distribution in a nontrivial way. This 
is an unavoidable numerical complication that we should face relentlessly if the analytical calculations are to be compared with 
Monte Carlo simulations. 

To this end, we have adapted a simulation method based on transition probabilities first proposed in P2J for SIR models 
in complex networks. The numerical approach considers all transitions between states that take place during the dynamical 
evolution of the subpopulations, defined by the system of differential equations ( |20] |. When dealing with structured populations, 
these transition rates depend, in general, on the connectivity class within which they occur Moreover, within each fc-class seven 
transitions are possible (see Fig.lTJ: 

- Birth of healthy individuals, 

- Natural death of healthy individuals, 

- Natural death of latently infected individuals, 

- Natural or disease-related death of sick individuals, 

- Transition from a healthy to the latent state, 

- Transition from a healthy to the sick (infectious) state, 

- Transition from a latent to the sick (infectious) state. 
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FIG. 2: Stationary proportion of sick individuals as a function of A (G [0.25,0.5]) for p = 0.07, r — 0.002, /itb ~ 0.2, /i — 0.009 and 
b = 0.01. The arrow marks the position of the epidemic threshold (Ac = 0.305) as given by analytical calculations using the previous values 
for the parameters and the degree distribution. Error bars are smaller than symbol size. The initial size of the network, characterized by a 
degree distribution P{k, to) ^ k^'^, is A^o = lO''. Each point is an average over at least 1000 values of (t)*. 



Each of these transitions is characterized by a characteristic transition rate cj^ ^ that can be directly derived from the system of 
equations that characterizes the rate at which they occur within the class k as; 

wi,fc = bNP{k,to) 
^2.k = liNP{k,t)uk 
uj3,k^ k^NP{k,t)lk 
W4,/c = (a* + fJ'tb)NP{k, t)tk 
i05,k = {l-p)XkNP{k,t)uke 
wg^fe = pXkNP{k,t)ukQ 
'ujT^k=rNP{k,t)lk 

Similarly, we define the sum of all these transition rates as the average rate at which one transition (of any kind) occurs: 

(24) 



This average transition rate in its turn defines the characteristic or average time r elapsed between any two consecutive transi- 
tions, the latter being defined as the inverse of 17: 

Given the previous definitions, the Monte Carlo algorithm is implemented in such a way that at each MC step (of duration r) 
one single transition takes place. Finally, the probability lii^k that a given transition actually happens, is calculated as: 

ni,/c = = TUJi,k, (26) 

that determines which of all possible transitions is realized at each time step r. 

We have made extensive numerical simulations of the model starting from an initial population made up of iV ~ 10^ in- 
dividuals, whose network of contacts follows an initial degree distribution P{k) ~ k~^. Moreover, every newborn joins the 
system with a degree that verifies the same connectivity distribution. As for the values of the parameters of the dynamics, and 
thinking of typical values for persistent diseases, we have set the following values: ji = 0.009 years"^, b = 0.010 years^^, 
pi'Tb = 0.200 years^^, r = 0.002 years"^, and p = 0.070. Demographical parameters b and fi are roughly those of a country 
like Spain, while the parameters p and r are in the range of typical values for the case of tuberculosis. has been chosen 
attending to numerical convenience (tuberculosis reaches a disease-related mortality rate that can be as large as 0.8). On the 
other hand, we note that a numerical criterion to define stationarity should also be adopted. In our simulations, we first let the 
system evolve for 4500 years and later take averages in a window of 10® Monte Carlo steps (which corresponds, roughly, with a 
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TABLE I: Dependency of the epidemic thresholds on the initial size of the network. Analytical values are obtained from Eq. l l23t using the 
moments of the distribution generated numerically. 





Analytical 


Numerical 






Aciv 


1000 


0.52903 


0.46968 


10000 


0.42559 


0.37595 


50000 


0.37620 


0.33088 


100000 


0.35854 


0.31926 



temporal lapse of 100 years), for the mean densities of healthy, latent and sick individuals defined in (|2TI) . This is a long enough 
time and ensures that all of the outputs to the left of the threshold are stabilized to the state (1, 0, 0) (this is achieved almost 
surely akeady for t < 4000 years). 

With the values for the parameters as specified above, the epidemic threshold is Ac = 0.305. In Fig. |2]i, we have plotted 
the stationary proportion of sick individuals for values of the probability of transmission A in the interval [0, 0.5]. As can be 
seen from the figure, the numerical and analytical results are in a reasonable agreement, despite the numerical challenges of 
simulating an open system in which the dynamics evolves on top of a complex topology. This indicates that the numerical 
method is accurate enough as to be used in situations where analytical predictions are not at hand. 

However, it is possible to carry out numerical simulations using a variation of the algorithm above in order to improve the 
accuracy in the determination of the epidemic threshold Ac. In this variation, instead of using a criterion for stationarity, we 
focus our attention in the vicinity of the critical value. More specifically, we first evaluate analytically the value for Ac and start 
the simulation there (obviously, if we don't have an analytical hint, the simulation can be started at any value of A). At each 
realization, we expect 6000 years for an eventual arrival of the system to the state {(u), (1), {t)) = (1, 0, 0). In the case that this 
state is not reached in that time, we assume that we are to the right of the critical point and so, we move to the left in A just a 
little quantity SX = 0.01. If, on the contrary, a minimum number of realizations (we used 10 in our simulations) stabilize at 
state {{u), (/), {t)) = (1, 0, 0), we assume to be to the left of the critical point and, consequently, we perform a d\ switch to the 
right. Each time that such kind of flip-flop algorithm (a sort of bisection method) changes direction, we divide by two the value 
of SX until the desired precision in Ac is obtained. Using this numerical approach, we have numerically calculated the values of 
Ac for different system sizes. The results are reported in Table |T] As expected from finite-size effects, the larger the size of the 
population, the smaller the absolute error between numerical and analytical thresholds is. Moreover, the larger the system size 
the smaller the epidemic threshold, which eventually should vanish in the thermodynamic limit. 



IV. CONCLUSIONS. 



We have discussed a model for the spreading of persistent infections in complex heterogeneous populations. The framework 
extends the epidemiological picture proposed in previous works. Our approach is particularly suited for diseases like Tubercu- 
losis, which shows large latency periods. The latent period results from the dynamical equilibrium that is established between 
the bacterium and the host's immune system, so that the host might not become infectious during its lifetime. These particular 
features makes it compulsory to work with an open system where newborns are continuously introduced in the population and 
individuals might die due to causes different from the disease itself. By assembling a model with all these ingredients, we have 
shown analytically that the epidemic threshold is proportional to the ratio between the first and second moments of the degree 
distribution. Therefore, our results point in the same direction that those obtained for the SIS and SIR model on top of the same 
topologies -the virtually unbounded connectivity fluctuations play a key role in the infection dynamics enhancing the epidemic 
incidence and lowering the epidemic threshold. 

From another point of view, we have developed a method well suited to numerically explore the dynamics of the system under 
study. In particular, we have been able to deal with the new challenge of having a system where the number of individuals in 
the population is not constant and, moreover, are connected following a given degree distribution (the well-mixed case is not 
challenging). Although we have applied the numerical approach to our particular system, it is worth stressing that it is general 
and can be applied to any problem for which transition rates between different classes and states are known. The results obtained 
agree well with the analytical estimates with the additional advantage that the whole phase diagram can be explored. 

Finally, it is also worth mentioning that the model discussed here is probably the simplest one may devise for the spreading of 
persistent infections in structured populations. However, despite the recent progresses in modeling disease contagion dynamics 
and pandemic outbreaks, the kind of spreading phenomena analyzed here is one of the issues that have remained less explored. 
Our aim is to take a first step towards more realistic modeling of persistent diseases. We have left for future investigation possible 
extensions of the current model that take into account the influence on the dynamics of vaccination, prophylaxis and recovery 
rates, as well as the effects of genetic heterogeneity in the pathogen but also in the host ll24ll . 
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